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ABSTRACT 



This thesis investigates the computer aided measurement of fiber 
diameters by laser diffraction. The proposed system consists of a light 
sensitive Random Access Memory (RAM) chip which collects light intensity 
data from the laser diffraction pattern. Measurements of the spatial 
location of the nodes of the diffraction pattern enables the calculation of 
the fiber diameter. These measurements may be performed manually which 
is tedious and requires subjective judgement of the nodes. The alternative 
method of direct processing of the intensity pattern was investigated. 
Simulation is conducted to examine the feasibility of this method. Results 
show such a system to be capable of providing one order of magnitude 
greater accuracy than optical microscopy measurements (with a shearing 
eyepiece) and double the accuracy of manual laser diffraction methods 
with the added advantage of permitting the option of total computer 
automation in data interpretation. 
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I. IMTRODUCTION 



Fiber reinforced composites are replacing structural metal 
components in today’s aircraft. The high strength and reduced weight of 
composites results in less drag, increased payload, and longer fatigue life. 
Additionally, the directional properties of composite materials provide 
unique design advantages over conventional materials. 

As with many developing technologies, the reasons for the success of 
fiber composites was not initially appreciated. Tsai [Ref. 1:p. 21 states 
that fortunately the modem composite was so strong it was reliable and 
competitive in spite of less than optimum design practice. Over the last 
twenty years, much progress has been made in understanding the micro and 
macro mechanics of composite materials. As this knowledge matures and 
is incorporated into the design process, the full potential of these 
materials can be realized. 

An important contribution to the design process is modeling 
structural reliability. Development of probabilistic models for a 
composite must take into account the complex relationships that exist 
between fiber and matrix. 
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One important parameter in the probabilistic model is fiber diameter, 
since variations in fiber diameter affect fiber failure density. Therefore, 
more accurate fiber diameter measurement results in an enhanced 
reliability model and a more quantitative prediction of structural 
reliability. 

The purpose of this research is to investigate a computer aided 
method of fiber diameter measurement. The existing procedure of 
diameter measurement by laser diffraction is accurate to within 0.5% 
[Ref. 2;p. 2101. It is desired to improve this accuracy by interpreting the 
diffraction pattern with a light sensitive RAM chip. 

The presentation begins with a discussion of diffraction pattern 
analysis. This is followed by an introduction to the Micron Eye’s theory 
and operation. Next, the simulation is discussed in depth followed by the 
results and conclusions. 
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II. BACKGROUND 



Fiber diameter measurement by laser diffraction is conducted by 
obstructing a collimated laser beam with a fiber sample. A diffraction 
pattern results and is characterized by alternating maxima and minima 
symmetric about a central maximum. (See Figure 1) 




Figure 1. Fiber Diameter Measurement by Laser Diffraction 



In previous work. Perry, Ineichen, and Eliasson, [Ref. 31 and Dennett, 
[Ref. 4], interpretation of the diffraction pattern consisted of finding the 
distance between interference nodes (minima). This distance is related to 
the fiber diameter. 
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Bennett used the slit approximation (Appendix A) to relate the 
distance between nodes to the fiber diameter. Perry, et a1., introduced a 
more exact solution by Kerker, [Ref 5:p. 260] and compared it to the slit 
approximation. Perry, et al., concluded that the slit approximation should 
be treated with caution. Therefore, Kerker's solution has been adopted for 
this study. 

In his paper, Bennett successfully demonstrated the feasibility of 
using the light sensitive RAM chip for diameter measurements (Ref 4]. An 
IS32 OPTICRAM MICRON EYE, manufactured by Micron Technology, Inc., 
was connected as a peripheral device to an Apple 11+ computer. By 
positioning the Micron Eye so that two interference nodes fit on its 
surface (see Figure 2), an "exposure" of the diffraction pattern could be 
printed by the computer (see Figure 3). Analysis of the printed diffraction 
pattern gave the diameter of the fiber using the slit approximation of 
Fraunhofer diffraction theory (see Appendix A). 




Figure 2. Micron Eye Array with Two Interference Nodes 
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Figure 3. Typical Micron Eye Photograph of Interference Nodes 
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III. DIFFRACTION PATTERN ANALYSIS 



Measuring fiber diameters by laser diffraction, requires an 
understanding of the diffraction pattern. This chapter introduces some 
features of diffraction patterns which will be useful in later analysis. 

A. FUNCTIONAL DESCRIPTION OF THE DIFFRACTION PATTERN 

The diffraction pattern of a fiber can be described as follows: 

l/l^= (2/K^LU) I ♦ 2 2 cos(ne) p (1) 

where 0 = the scattering angle 

= 2if /A ( A = laser wavelength) 

b„= J,(») / H,(2 )(q,) 

and oc- nd^A (d^ = fiber diameter) 

J^(ct) are Bessel functions of the first kind, 

H (2)(^) are Hankel functions of the second kind. 

II 
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A formal introduction of equation (1) and the related Fraunhofer 
diffraction theory is presented in Appendix A. 

Figure 4 depicts the three dimensional nature of the diffraction 
pattern. Equation (1) describes the intensities along one linear position, 
or slice, of the three dimension pattern. 



I 
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Figure 4. Three Dimensional Diffraction Pattern 



A plot of equation (1) is the Intensity Profile which shows the 
Intensity Ratio (1/1^) versus the angle theta, for a given fiber diameter 

and a given fiber to screen distance L. Figure 5 shows a typical Intensity 
Profile for a fiber diameter of 6.5 microns. 
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B. FEATURES OF THE INTENSITY PROFILE CURVES 



The intensity profile has several features worth noting. The first is 
the central (or zeroth order) maximum. The central maximum is the 
largest peak in Figure 5 and it is many times more intense than the 
subsequent maxima. It should also be obvious that the profile is 
symmetric about the central maximum. 

Other features of the curves are shown in Figure 6. Figure 6 is an 
enlarged view of one of the higher order maxima and associated minima. 
As discussed in Appendix A, the higher order maxima are not centrally 
located between the minima, rather they are displaced slightly towards 
the central maximum. Therefore, the higher order peaks are asymmetric 
about their maxima. This asymmetry means the maximum derivative on 
the upslope side of the curve will be greater than the maximum derivative 
on the downslope side of the curve. This difference is dealt with later 
when tuning for the optimum exposure is discussed. 

The minima in Figure 6 are also distinctive features of the curves. In 
Appendix A, it is shown that the location of the minima are explicitly 
related to the diameter of the fiber using the slit approximation; 

m = * of interference node 

®min ~ ^ mA./d 1 X - laser wavelength (2) 

d = diameter of fiber 
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Figure 6. Features of the Intensity Profile Curve 
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Equation (2) is only an approximation for a fiber, but it is useful for 
preliminary analysis. Such analysis includes: 

(1) determining where to place the Micron Eye for data collection and, 

(2) providing and initial guess of fiber diameter from diffraction 
pattern data. 

C. EFFECT OF DIAMETER ON THE INTENSITY PROFILE CURVES 

Figure 7 shows the effect of diameter on the Intensity Profiles. Note 
that successive maxima and minima are further from the central maximum 
for fibers of smaller diameter. This behavior will be important to later 
analysis. 

D. THE MAXIMUM DERIVATIVE AND THE MINIMUM NUMBER OF POINTS 

The goal of automated fiber diameter measurement is to find the 
diameter rapidly and accurately. Speed and accuracy can be exclusive. One 
can imagine that collecting all the points on the Intensity Profile cune 
will result in an almost exact determination of fiber diameter, at the 
expense of time. Two questions must be answered. First, what is the 
minimum number of points to uniquely describe an intensity profile? 
Second, where on the curve should these points be obtained? 
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Figure 7. Effect of Diameter on Intensity Profiles (3-D) 
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1. Minimum Number of Points 



Equation (1) consists of an infinite series. Determining the 
minimum number of points to uniquely define an infinite series is not a 
straightforward procedure. Cursory inspection of Figure 8 leads to the 
conclusion that one point will not uniquely define the curve, since 
intensity curves for many diameters pass through the same point (e.g., at 0 
= .04 radians). If the absolute intensity is not known, 2 points will not 
provide a unique solution either. For the case of non-absolute intensity 
measurements (as is the case in this experimental set-up) 3 points 
appears to uniquely define a curve. 

The curve in Figure 9 illustrates this observation. An Imaginary 
line Is fixed through three points (i.e., the experimental measurements). 
The absolute intensities of these three points are not known, and neither 
are the absolute spatial locations with respect to the central maximum 
(i.e., 6). Only the relative spatial locations among these points are 
known. Iteration of the diameter is equivalent to moving this line (with 
the x marks fixed relative to each other) along different locations on the 
curve, trying to match all three points. This iteration can be repeated on 
curves of different diameters (Figures 10 and 1 1) and no other match can 
be found. Thus, this argument defines three points as the minimum 
required to uniquely determine the intensity profile curve. 
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Figure 8. Effect of Diameter on Intensity Profiles (2-D) 
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Figures 9, 10, 11. Three Points to Define a Curve 
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2. Maximum derivative 



Having heuristlcally established that three points are required, the 
next step is to determine the optimum location from which to select the 
points. For example, if points are selected at the mimima (or maxima) a 
wide variation in theta results from a small change in intensity Ratio. 
This variation is defined by the derivative of the Intensity Ratio with 
respect to Theta, and at the extrema this derivative Is very small. It can 
be shown that the points on the curve where the derivative Is a maximum 
will have the least variation error. Therefore, it is desirable to use the 
Intensity Ratios corresponding to the maximum derivatives as the 
optimum Intensity ratios for the three points. These Intensity Ratios will 
be referred to as the Threshold Intensity Ratios. 

Figure 12 shows a portion of an Intensity Profile with a derivative 
curve (absolute value) superimposed. Vertical lines drawn through the 
derivative curve maxima intersect the Intensity Profile showing the 
location of the optimum (threshold) intensity points. 
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Figure 12. Intensity Profile with Superimposed Derivative Curve 
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IV. MICRON EYE 



The MICRON EYE image sensor Is an optically sensitive Random Access 
Memory (RAM) chip capable of sensing an image and translating it to 
digital computer compatible signals. The Micron Eye was selected to 
introduce automation to the existing techniques of fiber diameter 
measurement by laser diffraction. Automation is expected to increase the 
speed and accuracy of diameter measurements. This chapter introduces 
the theory and operation of the Micron Eye. 

A. PHYSICAL LAYOUT AND DIMENSIONS 

The Micron Eye (IS32 OpticRAM has two arrays each containing 128 
rows X 256 columns of sensors. This application will use only one of the 
arrays. The size of an array is 4420 jim by 877nm. (See Figure 13) 



4420um 



877|lm 

128rovs 




Figure 13. Micron Eye Array 
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Each sensor is a light sensitive element is called a pixel. The 
physical organization of the pixels is shown in Figure 14. The 128 x 256 
elements actually map into a 129 x 514 "cell placement grid". This 
arrangement leaves "space pixels" in between each pixel in the rov/ 
direction. The space pixels can be set high, low, or to the level which 
agrees with the majority of its nearest neighbors. 

B. THEORY OF OPERATION 

The pixels are capacitors which discharge a preapplied voltage at a 
rate proportional to both the intensity and duration of the impinging light. 
The voltage in an exposed capacitor is read and digitally compared to the 
fixed threshold voltage. If the voltage is below threshold the pixel is read 
as WHITE. If the voltage is above threshold the pixel is read as BLACK. 
The digital comparison concept will be an important part of the simulation 
in the next chapter. 

After a pixel is read, the row containing that pixel is refreshed. 
Refreshing sets all pixels which are below threshold to 0 volts, and all 
pixels which are above threshold to +5 volts. [Ref. 6] 
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IS32 OpticRAM™ topological information 




Figure 14. Micron Ege Physical Organization 
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1. Operation of the Micron Eue 

Operation of the Micron Eye is simple. The current research 
configuration uses an Apple Macintosh (512K) computer. Control software 
is provided by Micron Technology, Inc. The Micron Eye is connected 
through either the modem or printer port. (See Figure 15) The computer 
also acts as the power source for the Micron Eye. 




Figure 15. Macintosh and Micron Eye 



2. Exposure 

An exposure of a portion of the diffraction pattern can be made by 
varying exposure times. A sample exposure is shown in Figure 16. It is 
important to visualize that this exposure represents a cross section of the 
intensity profile curve. 
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Figure 16. Exposure and Its Relation to Intensity Profile 




Longer and shorter exposure times will vary the size of the cross 
section. Varying exposure time is equivalent to moving up or down the 
intensity ratio axis of the intensity profile curve. 

C. INTENSITY MEASUREMENTS WITH THE MICRON EYE 

Bennett [Ref. 4] used the Micron Eye to measure the distance between 
interference nodes. This approach was compatible with the slit 
approximation requiring only the theta locations of the interference nodes 
to give the diameter. In Kerker’s equation [Ref. 5;p. 260), one cannot solve 
explicitly for d. The solution requires knowledge of the Intensity Ratios 
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and their respective theta locations. For the Micron Eye, Intensity Ratios 
correspond to the user controlled exposure time. 

1. Two A p proaches to Intensity Calibration 

The Micron Eye must be calibrated to a reference intensity level. 
One method to accomplish this would require two exposures. One at + 

A and another at Ij ♦ A. Delta represents an unknown level above the 

zero (or absolute) intensity. (See Figure 17) The difference between 
these two values eliminates delta. 

A second method requires a fiber of known diameter. One exposure of 
this "calibration" fiber will allow determination of the absolute intensity. 



INTENSITY PROFILE 




Figure 17. Calibrating the Micron Eye 
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V. SIMULATION 



This chapter introduces the simulation of the Micron Eye / Laser 
Diffraction diameter measurement system. The simulation answers many 
questions relating to the real system. 

(1) where should the Micron Eye be placed? 

(2) what exposure is best? 

(3) what accuracy can be expected? 

(4) is the computer code valid? 

The simulation combines diffraction pattern analysis with Micron Eye 
operation. The worth of the results will depend on the accuracy of the 
simulation. 



A. COMPUTER PROGRAMS 

Several computer programs have been written to conduct this 
simulation. These are: 

(1) DATAMAKR 

(2) EXPOSURE 

(3) DIAFIND 



Appendix B discusses these programs in detail. 
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Briefly: 

(1) DATAMAKR is used to produce Intensity Ratio versus Theta data 
for a given fiber diameter d, and screen to fiber distance L. The 
user controls the number of points produce and their spacing. 

(2) EXPOSURE reads the Intensity Profile data generated by 
DATAMAKR and recommends the optimum exposure based on the 
average of the intensity ratios corresponding to the maximum 
derivatives. Exposure will also introduce random error into the 
data, as desired. 

(3) DIAFIND uses the data generated by EXPOSURE to find the diameter 
of the fiber. The programs begins with a user input guess of the 
diameter and conducts an iterative search / comparison routine to 
find the actual diameter of the data. 



B. OVERVIEW 

The general approach in examining the proposed system will be to consider 
a typical carbon fiber. This fiber has been assigned the arbitrary diameter 
of 7.254pim. Two additional fibers ±20% of the 7.254pim fiber are also 
considered. These fibers define the range of diameters from 5.803jjm to 
8.705pm. See Figures 18, 19, and 20 for Intensity Profiles curves for 
these diameters. 
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SIMULATION INTENSITY PROFILE 




Figure 18. Simulation Intensity Profile (5.803pm) 
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Figure 20. Simulation Intensity Profile (8.705)im) 



C. POSITIONING OF THE MICRON EYE 

In order to measure the diameter of the typical fiber, a method of 
postitioning the Micron Eye must be defined. Further, from this same 
position, it is desired to also measure the ±2035 fibers. This requires that 
one position of the Eye permit fiber diameter measurements over the 
specified range of fiber diameters. 
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Such a position is defined such that three data points will be obtained 
for any fiber in the diameter range, as in Figure 21. First, a simpler case 
will be discussed to introduce positioning for a given diameter fiber. 



Intensity Profile 



4 

at 



e 





3 Theta (radians) 



Micron Eye 

Figure 21. Exposure for Three Data Points 
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1. Positioning for a Single Fiber 

Suppose an exposure of the first through the second interference 
nodes (2”^ maximum) was desired for the 7.254iim fiber. Figure 22 shows 
the relationship between the fiber and the Micron Eye: 




Figure 22. Micron Eye Placement In Diffraction Pattern 



0^ = tan“* I X- / L 1 


(3) 


0H1 = t“"'' • *1.1 > L 1 


(4) 
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or tan 8i+i-tan 0i = 1/L (xi+i - xi) = Ax/L (5) 

lake AX = 4.4 mm (width of Micron Eye array) 

With Ax fixed, one can adjust L so that any desired Shi ■ ®i 

fit on the Micron Eye's array. Recall that = sin I mA/ d 1 

equation (2), which will give the theta locations of the first and second 
interference nodes. 

0 . = sin (mA/d) where m = node * 

mill 

A = laser wavelength (632.8nm) 
d = fiber diameter 

for 7.254jim; 0. = .08712 radians 

0|^j = . I 7459 radians (for i : I) 

with 0- and 0.^^ fixed, L is defined from equation (5); 

L = Ax / I tan 0.^1 - tan 0- ] = 49 mm 
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L is the distance of the Micron Eye fronn the fiber in the 
longitudinal direction. Referring back to equation (3), the lateral 
placement of the Eye is given (see Figure 23); 

^nner = ^ tan(0) = 4.27 mm 




y 

''inner 

I ► 

^otfter 

Figure 23. Micron Eye Position 



Thus, equations (3) through (5) give the Micron Eye position for any 
desired d, 0., L and m. 
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This position is not yet optimum for single diameter case. This 
analysis has captured only two interference nodes and the intervening 
maximum. Any exposure will yield information somewhere between the 
two extremes shown in Figure 24. The Micron Eye images correspond to 
the shaded areas of the intensity profile curve on the nght side of Figure 
24. 





Figure 24. Short and Long Exposures 



One method of finding the points would be to search the Micron Eye 
array data for a major axis. The endpoints of this major axis are the 
points required for the analysis. (See Figure 25) 
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Major Axis 




Figure 25. Endpoints of Major Axis 



Recall from diffraction analysis (Chapter 3) that three points are 
required to define the intensity profile curve. Therefore, more of the 
diffraction pattern must be seen by the Micron Eye. A reduction in L will 
yield the third point. (See Figure 26) L should not be reduced any more 
than is necessary, since fewer pixels are being used to describe the data, 
degrading the resolution of the Micron Eye. 




^ Major Axis 



endpoints 



Figure 26. Three Points on Major Axis 
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2. Positionino the Micron Eue for a Range of Diameters 

Consider the positioning of the Micron Eye so that three data points 
can be collected for a range of diameters. A "window" for the array must 
be defined which will ensure three points of data for any diameter in the 
specified range. 

Examine the Relative Intensity Profile plot in Figure 27. 0 begins 

at .075 radians which excludes the the central maximum. The central 
maximum is so Intense, the Micron Eye's fastest exposure cannnot prevent 
overexposure. Therefore, the central maximum Is not considered as a 
location for the window. 

The window should also be located in a region where the maximum 
intensity corresponding the the largest diameter fiber is close to the 
maximum intensity of the smallest diameter fiber. This ensures the 
Threshold Intensity Ratio will be more representative for all diameter in 
the range. This condition occurs at the higher order nodes (2 or greater), 
as illustrated in Figure 28. The plot of intensity ratio derivatives also 
shows this condition in Figure 29. 
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RELATIVE INTENSITY PROFILES 







Figure 27. Relative Intensity Profiles 
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I 




Figure 28. Region of Similar Maximum Intensities 
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To obtain three points, at about 1.5 maxima are required. Since 
the distance between nodes is the greatest for the smallest diameter, find 
the 6 locations for the smallest diameter which gives 1.5 maxima. 

sin (21/5.803|iin) = 216r>ds 
sin (31/5.803|un) = .321 rads 

Half of this range is about .050 radians and 1.5 maxima can be 

approximated by .165 > .321 radians. Admittedly, this is not a 

sophisticated method. The method is somewhat liberal and it was later 
discovered that a range of .165 to .300 radians was better than the larger 
range. This is because the minima are not of interest so the right hand 
minimum at .321 radians was discarded, and the range reduced to .300 
radians. Only the part of the curve where the slope is a maximum needs to 
be 'seen" by the array. 

A quick check can be made to insure that the cun/e representing 
the maximum diameter in the range fits inside the window, 
for 8.705pm; 

®min2 = sin (2 1 / 8.705) = . 1 50 rads 
^min3 = Sin (3 1 / 8.705) = .216 rads 
= sin (4 1 / 8.705) = .286 rads 



for 5.803pm: 

6 . ^ = 

min2 

0 - , = 
mm3 
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The next step is to determine the distance L at which the Micron 
Eye will cover the 0 range of .165 — > .300 radians. Equation (5) 

defines the longitudinal position and equation (3) defines the 
corresponding lateral position: 

L = Ax / llan(0-^^) - tan(0.)l = 30 mm 

X| = L tan(0|) s; 5 mm (inboard position of array) 
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VI. SIMULATION: THE PERFECT AND IMPERFECT DATA SETS 



Now that a window has been defined (0 range) for the diameter range, 
it is necessary to construct perfect data for the simulation diameters: d^, 
and dj. DATAMAKR generates this data for a given d and L, over the 
specified range of 0. Additionally, the number of points over the specified 
range must be chosen. The number of points defines the 0 interval and 
corresponds to the physical spacing of the pixels in the Micron Eye array 
which is on the order of lOjjm. To match this spacing, 400 points were 
generated over the interval .165 to .300 radians. 

A. DIGITIZATION OF DATA 

The EXPOSURE program is used to digitize the perfect data. Digitizing 
the data is the simulative analog of the Micron Eye threshold voltage 
comparison. In the simulation, any intensity ratios above a certain ratio 
(call it the threshold ratio) are assigned a value of 1 and any intensity 
ratio below the threshold ratio are assigned a value of 0. 
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1. Determining the Threshold Intensitu Ratio 

Before EXPOSURE can digitize the perfect data, the threshold ratio 
must be determined. EXPOSURE calls the subroutine DERIV which 
calculates the derivatives of the intensity profile data at each theta 
location. 

DERIV proceeds to search for all the local maximum derivatives. 
The average of the maximum derivatives is returned to the EXPOSURE 
program and is used as the threshold intensity. 

After the data is digitized with respect to the threshold intensity, 
the digital data is searched to find the theta locations where the intensity 
changes from 0 1 or 1 — > 0. The values for theta at these locations 

are averaged and this average is taken to be the theta location where the 
threshold intensity is located. The averaging of theta is based on the fact 
the Intensity Profile curve is approximately linear in the region of the 
Threshold Intensity Ratio. Figure 30 shows the general location of the 
Threshold Intensity Ratio. 

In the simulation, EXPOSURE finds three average theta locations 
corresponding to the threshold intensity ratio. These data are stored In a 
data file which is subsequently input to the program DIAFIND. DIAFIND 
prompts the user to input an initial guess of the diameter and proceeds to 
find the diameter associated with the EXPOSURE data. 
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B. INTRODUCTION OF ERROR INTO THE PERFECT DATA 

EXPOSURE prompts the user for an error input. If error is desired, 
EXPOSURE calls the subroutine RANDOM. RANDOM introduces error based on 
the maximum intensity ratio in the perfect data; 

U = 'e * I '„»c * error • RND 1 

Where ERROR is the error to be introduced and RND is a random number 
between +1 generated by the NONIMSL subroutine RANDU. RANDOM returns 
the now' "imperfect” intensity ratios to the EXPOSURE program. 

C. DIGITIZATION PROBLEMS WITH THE IMPERFECT DATA 

Digitizing the imperfect data results in local regions where the 
intensity oscillates between 0 and 1, as shown in Table 1. The 
introduction of random error changes the smooth curve to erratic points. 
Locally, variations above and below the threshold intensity ratio cause a 
series of 0-->1 and 1—>0 oscillations. (See Figure 31) Ideally, only one 
local theta location is to be associated with the threshold intensity ratio. 

There are two approaches to this problem. First, the average of the 
local theta locations can be calculated. This method results in three 
theta locations for input into the DIAFIND program (the same as perfect 
data). 
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TABLE 1. DIGITIZATION ERROR 



THETA LOCATION INTENSITY RATIO 



0. 2321624000E 00 
0. 2324997000E 00 
0. 2328373000E 00 
0. 2331748000E 00 
0. 2335124000E 00 
0. 2338498000E 00 
0. 2341873000E 00 
0. 2345249000E 00 
0. 2348623000E 00 
0. 2351998000E 00 
0. 2355374000E 00 
0. 2358747000E 00 
0. 2362123000E 00 
0. 2365499000E 00 
0. 2368872000E 00 
0. 2372248000E 00 
0. 2375622000E 00 
0. 2378997000E 00 
0. 2382373000E 00 
0. 2385748000E 00 
0. 2389124000E 00 
0. 2392498000E 00 
0. 2395873000E 00 
0. 2399249000E 00 
0. 2402623000E 00 
0. 2405998000E 00 
0. 2409374000E 00 
0. 2412747000E 00 
0. 2416123000E 00 
0. 2419497000E 00 



0. 3544308000E-04 1 

0. 3196243000E-04 1 

0. 3160309000E-04 1 

0. 3246925000E-04 1 

0. 3372986000E-04 1 

0. 3785736000E-04 0 

0. 3090432000E-04 1 

0. 3342219000E-04 1 

0. 3387519000E-04 1 

0, 3205373000E-04 1 

0. 3516222000E-04 1 

0. 3860490000E-04 0 

0. 3421140000E-04 1 

0. 3468163000E-04 1 

0. 3580747000E-04 1 

0. 3904779000E-04 0 

0. 3578593000E-04 1 

0. 3912353000E-04 0 

0, 3599435000E-04 1 

0. 3474336000E-04 1 

0. 3538676000E-04 1 

0. 3470230000E-04 1 

0. 3952176000E-04 0 

0. 3923657000E-04 0 

0. 3944175000E-04 0 

0. 3966210000E-04 0 

0. 4409726000E-04 0 

0. 4269459000E-04 0 

0. 4371968000E-04 0 

0. 3904802000E-04 0 
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INTENSITY RATIO *10 

27.0 29.0 31.0 33.0 35.0 



ERROR SIMULATION 
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The second method is to accept the theta locations as they are and 
enter them into the DIAFIND program. The latter method will pass a 
variable number of data points to DIAFIND. 

Both approaches were tested (at 5^ error for 7.254 pim) and the same 
diameter was recovered by DIAFIND in each case. The only difference 
between the two methods is that the averaging method runs one second 
faster (out of 22 seconds on the IBM 360) than the other method. 

The averaging method was adopted for this simulation for two 
reasons: 

(1) The data input to DIAFIND will always be a constant number of 
points so that DIAFIND can be easily adapted to run actual Micron 
Eye data without simulation related logic buried in the code. 

(2) The averaging method returns the same diameter as the other 
method. 

D. CHOOSING THE CORRECT EXPOSURE (TUNING) 

The optimum exposure for this simulation was calculated by the 
program EXPOSURE. It was an average value of the intensity ratios 
associated with the maximum derivatives locations for the 7.254pim fiber 
data. This same exposure was used throughout the simulation for all the 
fibers. 

Using the same exposure introduces realism into the simulation since 
one cannot tune the exposure for eveny fiber that is to be measured. 
Because the Micron Eye window is located at least two nodes out from the 
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central maximum the exposure is closer to optimum for all diameters over 
the range. (In contrast to a window location closer to the central 
maximum.) 

The process EXPOSURE uses to determine the threshold intensity ratio 
can be called “tuning". The Micron Eye analog of tuning would consist of 
three steps; 

(1) Determine the approximate diameter for the fiber to be measured 
using an optical shearing eyepiece. 

(2) Determine the Micron Eye location parameters L and x by the 
methods in Chapter 5. 

(3) Vary exposure to obtain a defined relationship between the three 
points on the Micron Eye photograph. 

As an example, return to the case where L = 30 mm and x. _ = 5mm. 
(x^^ is the inboard x location of the Micron Eye array.) Figure 32 shows 

the "tuned" relationship between three points on the Micron Eye array. 
This relationship can be expressed as a ratio of a . b : c. The theta 
locations at the inboard and outboard edges of the Micron Eye are known, 
e g-, 0.,^ = .165 radians and 0„.. = .300 radians. 

^ * Tnn«^ outer 
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e. 

inner 



e, 



e. 



0 outer 



Figure 32, Three Point Spacing Ratio 



inner < 6| < 02 < ©2 ^ ®ouier * 

To find 6 ^ , Q^, and 63 . njn EXPOSURE for the diameter being tuned. 

EXPOSURE outputs the theta values for the optimum exposure. (See Figure 
33) 



The theta values are related to x locations by; 

X, = L tan(e,) 

X2 = L tantGj) 

X3 = L tanCSj) 

where a = x, - 

b : Xj -X, 

C = X3-X2 
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Figure 33. Finding the Theta Locations 



VII. DISCUSSION OF RESULTS 



The simulation was conducted for three diameters; 
d^ = 5.803 urn 

62 = 7.254 )im 

dj = 8.705 (where d^ and d^ are ±20% 62 ) 



For each three levels of error were introduced: 1%, 2% and 5%. 
Figure 34 shows the three intensity profile curves for the d^. 
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The simulated Micron Eye window was determined as previously 
described at L = 30mm and = 5mm. The exposure was tuned with 

respect to the 7.254pim perfect data. EXPOSURE recommended the 
Threshold Intensity Ratio of .0000370. This Threshold Intensity Ratio 
remained constant throughout the simulation, for all diameters. 

Before the simulation, all diameters were tested with no error using 
the 7.254pim Threshold Intensity Ratio. DIAFIND recovered d^, 62 and dj 

exactly (i.e., to the three digits accuracy of the original diameters). 

The simulation originally began by collecting thirty data points for 
each diameter/error combination. Twenty additional points were 
collected (50 total) to produce meaningful histograms. 

Figures 35 through 43 are histograms depicting the results of the 
simulation. In general, the expectation was to see a decrease in 
resolution as more error was introduced. It was also anticipated the 
resolution would be best for the 7.254pim case (for all values of error) 
since the exposure was "tuned” for this diameter. Finally, it was hoped the 
method would be more accurate than existing laser diffraction 
measurement methods (s 0.5%). 
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5.803 Jim fiber (1% error) 
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Figure 35. 5.803jim Fiber (IX error) 
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7.2461 



7.254 |im fiber (IX error) 




Diameter (microns) 

Interval = .0042tim 



Figure 36. 7.254jim Fiber (IX error) 



7.2629 
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8.705 Jim fiber (IX error) 




Diameter (microns) 

interval = .006jim 



Figure 37. 8.705jim Fiber (IX error) 
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5.7884 



5.803 urn fiber (2% error) 




Diameter (microns) 

interval .0048um 



Figure 38. 5.803um Fiber ( 2 % error) 
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7.254 |im fiber (2X error) 
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Figure 39. 7.254jim Fiber (2% error) 
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Figure 40. 8.705|im Fiber (2X error) 
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5.803 nm fiber (5% error) 
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Figure 41. 5.803jim Fiber (5X error) 
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7.254 iim fiber (5% error) 



22 

20 

18 

16 

^ 14 

S 12 

3 

9 10 

O 

^ 8 

6 

4 

2 

0 



m 


m 


iO 


10 


to 


10 




vD 


K) 


O 


h- 




N 


N) 




n 


lO 






04 























Diameter (microns) 

interval = .007)jm 



Figure 42. 7.254|im Fiber (5% error) 
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7.2715 





8.705 ym fiber (5X error) 
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Figure 43. 8.705pm Fiber (5% error) 
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A. ACCURACY VERSUS RESOLUTION 



The histrograms show the accuracy and the resolution of the method. 
The accuracy is associated with the largest spike, and related to the size 
of the interval called the resolution. The resolution is the ability of the 
routine to “see” a difference between two different fibers. For example, if 
the resolution is .OOSpim the method does not discriminate between 
8.705|im and 8.700pim. Table 2 shows accuracy versus resolution for the 
results. 



Table 2. ACCURACY VERSUS RESOLUTION 



TERROR 


DIAtlETERpim 


RESOLUTION 


ACCURACY (Res/dact) % 


1 


5.803 


.004 


.069 




7.254 . 


.0042 


.058 




8.705 


.0060 


.069 


2 


5.803 


.0048 


.083 




7.254 


.0054 


.074 




8.705 


.0049 


.056 


5 


5.803 


.0105 


.181 




7.254 


.0070 


.096 




8.705 


.0136 


.156 
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The data show a decrease in accuracy and resolution as error is 
Increased. It is evident in most cases that the 7.254|im results are 
better because it was the “tuned" diameter. The only exception is the 
accuracy for the 8.705tim fiber ( 2 % error) is better than the 7.254jim 
fiber. Overall, the largest error is .18 percent which is less than one half 
of the error associated with the manual laser diffraction method. 

8. ERROR 

There are two contributions to error in this simulation. 

The first can result from programming/calculation errors. Many trial 
runs of the software were taken to minimize the likelihood of this kind of 
error. 

The next type of error is the digitizing error encountered in an actual 
experiment. This error is simulated based on the maximum intensity of 
the profile curve, in the region of interest, to apply random error equally 
to all points. Had the random error been based on each point, points with 
less intensity would have less error, and points with higher intensity more 
error, which is inconsistent with the physical digitizing process. Using 
the maximum intensity avoided such a condition but it cannot be 
ascertained that this is the optimal representation of the physical system. 
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VIII. CONCLUSIONS 



The results of the simulation are encouraging. The largest error is 
.18^ or less than one half of the manual methods. Because the programs 
were carefully developed and tested it is unlikely they contributed to the 
error. Also, much effort was directed towards accurate simulation of the 
physical system so that the results would reflect what can be expected 
from the actual experimental measurements. 

The simulation demonstrates an increase in accuracy two to ten 
times better than that currently possible by making manual measurements 
with laser diffraction. The method also lends itself to automation which 
makes it attractive for quality control purposes and research for 
materials development. 
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(}(■ RECOnriENDATIONS 



This study has shown the feasibilty of computer aided diameter 
measurements. There are many directions future work can take. 

A. SOFTWARE. 

More development and testing of the software could result in 
increased accuracy. It would also be valuable to implement the softv/are 
on a small computer (like the Macintosh) to allow real time processing of 
actual data. 

B. HARDWARE. 

% 

There remain some hardware considerations which must be resolved. 
The biggest of these, perhaps, is the accurate positioning of the Micron Eye 
array in the direction perpendicular to the laser beam. In an effort to 
increase the resolution and accuracy of the problem, the system may 
benefit from two Micron Eyes. 
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APPENDIX A. FRAUNHOFER DIFFRACTION THEORY 



The following Is a brief discussion of Fraunhofer diffraction theory 
with emphasis on aspects of the theory which relate to this research. 

A. FRAUNHOFER DIFFRACTION 

Fraunhofer diffraction (Figure 44) results when light approaches and 
leaves an obstacle or aperture in the form of plane wavefronts. [Ref. 4:p. 
176] The light source and the plane of observation in effect are at 
infinity. A collimated laser beam is ideally capable of Fraunhofer 
diffraction because its beam consists of parallel rays advancing in phase. 



Central Maximum 



biterferenee Ntxles 




Figure 44. Typical Single Slit Diffraction Pattern 
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B. THE CLASSICAL SINGLE SLIT EXPERIMENT 

The simplest demonstration of Fraunhofer diffraction is the single 
slit experiment. (See Figure 45) Parallel, collimated light passes through 
a slit of width a. The diffraction pattern is visible on a screen located a 
distance L from the slit. An observer at point P, moving across the 
screen, sees a succession of maximum and minimum intensity points. 
These extrema are the result of constructive and destructive interference 
of the light. For example, a minimum occurs when the angle 0 produces a 
phase difference of one wavelength between the rays at the upper and 
lower edges of the slit. Thus, minima occur whenever 

a sin(0) = mA, where m = 1 , 2, 3 (A. 1 ) 

These minima are referred to as interference nodes. (For further 
discussion of this subject, see Meyer-Arendt, (Ref. 41.) 
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Figure 45. Single Slit Path Difference Relation 
to the Interference Node 



1. Diffraction Minima 

Examination of equation (A.l) shows that as the slit becomes 
narrower the angle 0 becomes larger. As this theory is extended to 
approximate the diffraction phenomena of an obstruction (a fiber), one can 
expect short distances between interference nodes for larger fibers and 
greater distances between nodes for fibers of smaller diameters. 



76 



Another important point concerns the distance between minima. It 
appears the interference nodes are equidistant. This is true only within 
the limits of the small angle approximation. The distances between 
interference nodes actually increases as one moves outward from the 
central maximum by the relation: 

^min = Sin I mA/d 1 (A.2) 

2. Diffraction Maxima 

The diffraction maxima are not located midway between minima. 
The locations of the maxima can be derived as by Meyer-Arendt [Ref. 7:p. 
220). These occur whenever the derivative of the intensity is equal to 
zero: 



dl/dB = dl/dB [ (sin B/B)^ 1 = 0 

= 2lo sinB/B I -sinB/B^ + cosB/B 1 

or whenever: B = tanB (A.3) 
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As per Houstoun (Ref. 8;p. 68) the maxima can be seen where 
y=tonB and y=B. (See Figure 46) The maxima are displaced slightly 
from center towards the central maximum. Much further away from the 
central maximum, the maxima are nearly halfway between minima. 




Figure 46. Locations of the Diffraction Maxima 
C. FRAUNHOFER DIFFRACTION AND THE FOURIER TRANSFORM 

The foregoing discussion centered on the Fraunhofer diffraction due 
to a slit. The mathematics of the slit example are simple and provide 
insight into diffraction physics. A more elegant derivation of Fraunhofer 
diffraction shows that the diffraction pattern of an object is the Fourier 
transform of that object. (Ref. 9 p. 174] 
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for the slit; 




Figure 47. Fourier Transform for the Single Slit 



Thus, F(0) describes the amplitude of the diffraction pattern, and 
IF{e)p represents the intensity. 

Modeling the transform for the obstruction is more complicated, 
consider the complement of the slit transform; g(x) = 1 - f(x). 




Figure 48. Complement of the Fourier Slit Transform 
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This results in the form I oo - 2sin (b6)/6 1 which is not 
transformable. 

An alternative is to use a substitute function h(x) = g(x) - f(x). 
Physically, this is approximating the obstruction as two parallel slits: 




Figure 49. Fourier Transform for Two Parallel Slits 

which has the solution: 

H(0) = (2sin c0)/0 - (2sin b0)/ 0 
where c would be large, but not infinite. 

This is still a crude approximation which shall be improved upon in 
the next section. 
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D. REJECTION OF THE SLIT APPROXIMATION 

In their paper titled Tiber Diameter Measurement by Laser 
Diffraction” [Ref. 3:p. 1378], Perry, Ineichen, and Eliasson conclude that 
the diffraction pattern of a real fiber is sufficiently different from that 
of a slit to warrant the slit approximation being treated with caution. 

Further, they recommend a solution presented by Kerker lRef.5:p. 260] 
which has been adopted in this study. Kerker’s solution assumes the fiber 
is perfectly reflecting (i.e.,the reflective index m=oo), and while this is 
not completely true. Perry, et al., {Ref. 3:p. 1378] indicate some degree of 
absorption is not likely to be significant. 

E. INTENSITY EQUATION FOR A REAL FIBER 

Kerker {Ref. 5:p. 260] gives the scattered intensity relation fora real 
fiber: 

I/I = (2/K Lif) I b ♦ 2 2 b cos(ne) P (A.4) 

O O On 

where 0 = the scattering angle 

= 2if /Tv ( A = laser wavelength) 

b^= J„(or) / 
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and 



o-r ud^A (df = fiber diameter) 

J (or) are Bessel functions of the first kind, 

are Hankel functions of the second kind. 

The real fiber equation (A.4) is somewhat obscure in its compact 
form. It can be shown that; 

l/l^=(2/K^LiJ) [( r^+22 r^cos(n0) )2+(s^+22 Sj^cosCnO) )^] (A.5) 

where ] (A.6) 

and 1 (A.7) 

Now, the intensity ratios for any 0 location can be calculated for any 
diameter fiber. 

1. Sensitivity of the Results to the Number of Bessel Terms 

The calculation of Intensity Ratios requires the computation of n 
and Bessel terms, equations (A.6) and (A.7). Here the phrase 

"Bessel terms" indicates the algebraic combinations of the and 
Bessel functions. 



82 



a. Two competing phenomena 

There exist two competing phenomena which govern the number 
of Bessel terms to be used in the calculations. The first requires a 
minimum number of terms for accuracy. The second, limits the number of 
terms so the Y functions do not cause an underflow error during 
computation. 

(1) Minimum Number of Terms . As in any series, there is a 
minimum number of terms required for computational accuracy. Figure 50 
shows the effect of the number of Bessel terms computed for. seven 
curves, all with a diameter of eight microns. 

Note that the 43, 50, 75, and 86 term curves are nearly identical and 
the curves with fewer than 43 terms are decreasing towards a smooth 
line. This research, indicates that 50 terms returns four digit accuracy 
for diameters ranging from 5 to 10 microns. 

(2) Maximum Number of Terms. The maximum number of terms 

depends on the value of oc. As the diameter decreases, so does or. This 
in turn produces very large values of (in the and 

denominator) which causes an underflow error. For the diameter range 
considered in this research, underflow was not a problem. In one 
instance, below 5 microns, it was found the number of Bessel terms had to 
be reduced to 43 to prevent underflow. The effects of the minimum 
number of terms was not investigated at this smaller diameter range. 
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INTENSITY PROFILE 
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Figure 50. Effect of Number of Bessel Terms 
2. Effect of Different Diameters 

Recall that the diffraction pattern for a slit showed that as the 
diameter decreased, the interference nodes moved away from the center of 
the pattern. It is interesting to note that the pattern for a real fiber 
behaves the same way. (See Figure 51) 
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INTENSITY PROFILES 




Figure 51. Effect of Different Diameters 
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APPENDIX B. COMPUTER PROGRAMS 



The following is a list of programs discussed in this Appendix; 

datamakr 

DIAFIND 

EXPOSURE 



The programs are written in Waterloo Fortran IV (WATFIV) and run on 
the Naval Postgraduate School's IBM 360 computer. 

A. FORMAT OF DATA 

All data is formatted using exponential notation. The most 
frequently manipulated files are those containing theta locations and 
intensity ratios. The format for these files is: (IX, El 7. 10, IX, El 7. 10). 
In WATFIV all data files must be of filetype “WATFIV". 

To compile and run a program on the IBM 360 type “ WATFIV 
PR06NAME DATAFILE *(XTVPE", where PR06NAME is the filename of the 
program and DATAFILE is the filename of the data file. 
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B. COMPUTER PROGRAMS 

This appendix presents the various computer programs and outlines 
their logic. 

1. DATAMAKR 

DATAMAKR generates the data used in the simulation. The equation 
for computing the Intensity Ratios of the diffraction pattern for any value 
of theta is; 

l/l^ = (2/K^Lfl) I + 2 5 cos(ne) P (A.4) 

where 0 = the scattering angle 

= 2if /A ( A = laser wavelength) 

and ckr= ud^/A (d^= fiber diameter) 

are Bessel functions of the first kind, 
are Hankel functions of the second kind. 

All computations begin by calculating the required number of 
Bessel terms. Two NONIMSL subroutines are called; BESJ for the J„(<ar) and 
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BESY for the Yj^Ca). Note: should double precision be desired, the NONIMSL 
subroutine BINT will return double precision values for the and Y^. 

The Bessel function values are stored in an vector J(l) and Y(l), 
and are used to calculate values of R(l) and S(l). Next, the matrix of 
cosines is contructed. The size of this matrix is M by K, where N = M-1 
and N is the number of Bessel terms. K is the number of theta locations. 
Values of theta at each location have previously been stored in an array 
called TOO. 

To visualize the program’s calculations, the matrices symbolic of 
these calculations are sketched. (See Figure 52) Note that and are 

outside the summation term in equation (A.l) and therefore will not be 
multiplied with any cosine terms. Since array subscripts must be 
designated with a nonzero value, and will be represented by r(l) 
and s( 1 ). 

After multiplying the two matrices, it is necessary to complete 
the summation by summing the columns in the matrix. This results in an N 
term value for each value of theta. A straightforward calculation then 
yields the intensity ratios. The intensity ratios are stored in an array 
with their corresponding theta locations. 
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Figure 52. DATAMAKR Matrix Algebra 
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The logic of DATAMAKR is; 



Input > * Bessel terms to be computed 

diameter of the fiber 
screen to fiber distance, L 
array of theta location values 

Calculate: and using BESJ and BESY 

Produce the Matrix of Cosines 
Produce the Matrix of Products 
Sum the columns of the Product Matrix 
Compute and output the Intensity Ratios 



2. DIAFIND 

This program finds the diameter of a fiber through an iterative 
process of residual comparison. The program accepts the data which has 
been output by the program EXPOSURE. The user is prompted for K and dj. 

K is the number of theta values the program should expect and d. is the 
initial guess of the fiber diameter (in microns). 
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The program first calculates intensity ratios at the diameter d., 
corresponding to the input theta locations. The difference between the 
input intensity ratios and those associated with dj is called the residual. 
The program calculates another set of intensity ratios based on a new 
diameter d, * Ad^, where Ad. is a small increment of diameter (usually 

.05 microns to start). The intensity ratios calculated using d. + Ad. are 

compared with the input intensity ratios to give a second residual. The 
residuals are compared and program logic determines whether or not the 
Ad increment is producing convergence to the actual diameter. The 
initial guess diameter is incremented and decremented as necessary until 
a desired level of accuracy is achieved. 

A key to understanding the convergence process is the residual curve 
in Figure 53. It is important that the initial guess, d., be fairly close to 
the actual diameter. 

For example, the residual curve in Figure 53 is for an actual diameter 
of 7 microns. If d. is greater than 9 microns, the program will not 

converge to the correct diameter, but rather to a diameter just above 10 
microns. 

It is hard to define exactly the limits of d- for any given actual 
diameter. Figure 54 is a residual cun/e for 5 microns and shov/s an upper 
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limit of 9 microns for dj. A thorough study to define limits for dj has 
not been conducted, although convergence has always been attained by 
guessing d. within +1 micron for diameters ranging from 5 to 9 microns. 



RESIDUAL VS DIAMETER 



a 

o = 




Figure 53. Residual Versus Diameter for 7|im Fiber 
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RESIDUAL *10 

10.0 15.0 20.0 25.0 30.0 



RESIDUAL VS DIAMETER 




Figure 54. Residual Versus Diameter for 5fim Fiber 
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3. EXPOSURE 



EXPOSURE simulates the operation of the Micron Eye by taking a 
photograph of intensity ratio data. EXPOSURE reads in the perfect data 
generated by DATAMAKR and digitizes the data vyith respect to a threshold 
intensity ratio. 

The optimum intensity ratio for an exposure is termed the 
threshold intensity ratio. The threshold intensity ratio is located 
where the absolute value of the derivative of the intensity ratio v/ith 
respect to theta is a maximum. This assures that the threshold intensity 
ratio is located in a region of the curve which is closest to a straight line. 
This reduces the effects of subsequent interpolation errors. 

The threshold intensity ratio is calculated by the subroutine 
DERIV. Since the input intensity profile curve will have at least three 
maximum derivative points, DERIV calculates the average of the three 
intensity ratios. This average is then passed to the calling program, 
EXPOSURE. 

EXPOSURE searches the intensity ratios and compares them with 
the value of the threshold intensity ratio. This process corresponds to the 
Micron Eye addressing a pixel and comparing its voltage to the threshold 
voltage. If an intensity ratio is greater than the threshold value, the ratio 
is assigned a digital value of 1. If the intensity ratio is less than the 
threshold value, the ratio is assigned a digital value of 0. 
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The digitized data is then searched to find the theta locations at 
which the digital intensities change from 0 --> 1, or from 1 — > 0. The 
theta locations are averaged to produce a theta location which is very 
close to the threshold intensity. Because this averaging process is 
occuring in the regions of steepest slope, error is assumed to be 
minimized since the curve can be approximated as a straight line (refer to 
Figure 30 in Chapter 6). 

EXPOSURE provides the user with the option of introducing error 
into the perfect data. The program prompts the user to input the desired 
error and calls the subroutine RANDOM. RANDOM introduces random error 
into the intensity ratios and returns the imperfect data to EXPOSURE. 
EXPOSURE digitizes the imperfect data and searches for the occurences 
where 0-->1 and 1-->0. Imperfect data will have many more occurences 
of the digital intensities changing from 1— >0 and 0-->1. EXPOSURE will 
average the theta locations if they are in the same location (i.e., v/ithin 
±.005 radians). This method has been tested for random error up to 5 
percent. 
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CALCULATE THE R(N) AND S(N) 
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PROMPT USER FOR INPUT INITIAL GUESS DIAMETER: 













lO 
























. 
























CM 


- 






















rH 
























u 


CO 










- 














p 










Eh 












- 


p 










CO 






hH 






cn 


EH 










p 




in 








HH 


w 










H 




Z 










S 














0 


< 






in 
















0^ 


EH 






P 


- 














0 


< 






w 
















HH 


Q 






Eh 


10 














S 








W 


. 






















S 


CM 














IS 


Eh 








rH 














hH 


HH 






p; 


w 










P 






3 : 






HH 












C 














- 




P 






EH 




pc^ 


P 






EH 












P 




U 


W 






p 


CO 










P 




H 


Eh 






p 


HH 




EH 






P 




W 


< 






S 
















S 


hH 






HH 


P 














< 


0 








0 




S 










HH 


0 






P 


S 














Q 


in 






W 


< 
















in 






Eh 


Eh 




X 










CO 


< 






U 


CO 










P 




in 






,, 


S 


HH 










P 




w 








< 


P 




in 










ID 


w 




w 


HH 






p 






P 




0 


H 




w 


P 


P 




0 






S 






w 




pci 




w 




Eh 






C 




j 






0 


CO 


p 




IS 






HH 




< 


< 




in 


CO 


HH 




HH 






P 




hH 


HH 






w 


P 














EH 


Q 




w 


p 






V 










HH 






K 


0 


0 




EH 










S 


P 




EH 




EH 




0 










HH 


< 






- p 






C 












D 




0 


< 


s 




EH 








2; 


EH 


EH 




EH 


S hh 


p 




S 






P 


pc^ 


ID 


0 






P--B 


p 




HH 






HH 


0 


Pk 


< 


VO 




OrH HH 


p 










CO 


W 


s 




0 


EH 


C 0 w^ 2 ; 


0 










P 




^ hH 


EH 


1 


D 


PUhh 


in 


VO 


P 






P 




- - 


P--- 


p 


P 


PW- 


- 


0 


P 








0 




P-r- 


iH 


S 


OP - 


P N 


1 












- X 


S . 


W» 


HH 


- X 


X 


p 










«rH 


— *-» ,P 


< 




1 


rH 


LOO 


P 








W(N - 


LD- .< 


^ HH 


0 


COvO ^ 


*- 


OCM 








HH 


SrH- iH- 


rH ^ .D 


P 


P 


SrH> 


rH- 


. 1 




0 




- * 


0 


iH \< 


X .EH 


< 


0 


0 \ 


\ 


P 


rHHHPt^ 


rH 


oocr> 






^ hH 




HH 11 


w 


Eh — 




0 11 rH 


wCOP 


II 


rH (Jl 


> — — 


OhVOwaDwO 


vDwUK 


P^ 




pvOw^ 


VO*w 


. 


PP 


II 


II II 




C14 '‘w^H "w^H 


^wfH2; 






P'wtH 


^wtH 


OP II 


PP II 


P 


11 II 


V ^ ^ 


WCClK - 


wco - 


11 rH 




P< 


PC 


c 


P 


PEH 


PP 


EhEhEh 


JHSHSQ 


EhSP:^Q 


— . 




PEhS 


EhS 


II HC 


P'^ 


PS 


OPPEH 


sss 


JhHpriHHpCiC 




<u 




pHHpc:J 


HHpci 


PEH 


II PrH 


SDrHCSCkW 


HHHHHH 


CpciOPciOW 


Pi^OHHU 


HHU 




<PO 


PO 


PPP 




CdO II iJOCMi] 


PPP 


o:^Ci4 2:Ci4P^ 


3:Cx4SPc^ 


PP 




OS:P 


3 :P 


PPP 


POP 




PPP 


























iH 


cs ^ 


LO 






VO 














iH 


rH rH 


r— 1 






1 — i 


rH 













r— i iH iH 

o o ooo oo oo ooo oo o 



106 



WHILE( ABS(TEST) . GT. BETA) DO 
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GENERATE THE MATRIX OF PRODUCTS 
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DO 900 1=1, K 

INTDIF(I) = ( INTGES(I) - INTACT( I ) 
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